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Abstract 

We have performed a numerical investigation of the ground state properties 
of the frustrated quantum Heisenberg antiferromagnet on the square lattice 
("Ji — J2 model"), using exact diagonalization of finite clusters with 16, 20, 
32, and 36 sites. Using a finite-size scaling analysis we obtain results for a 
number of physical properties: magnetic order parameters, ground state en- 
ergy, and magnetic susceptibility (at g = 0). For the unfrustrated case these 
results agree with series expansions and quantum Monte Carlo calculations 
to within a percent or better. In order to assess the reliability of our calcu- 
lations, we also investigate regions of parameter space with well-established 
magnetic order, in particular the non-frustrated case J2 < 0. We find that 
in many cases, in particular for the intermediate region 0.3 < J2/J1 < 0.7, 
the 16 site cluster shows anomalous finite size effects. Omitting this clus- 
ter from the analysis, our principal result is that there is Neel type order 
for J2/J1 < 0.34 and collinear magnetic order (wavevector Q = (0,7r)) for 
J2/J1 > 0.68. There thus is a region in parameter space without any form of 
magnetic order. Including the 16 site cluster, or analyzing the independently 
calculated magnetic susceptibility we arrive at the same conclusion, but with 
modified values for the range of existence of the nonmagnetic region. We also 
find numerical values for the spin-wave velocity and the spin stiffness. The 



spin-wave velocity remains finite at the magnetic-nonmagnetic transition, as 

expected from the nonhnear sigma model analogy. 
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I. INTRODUCTION 



In this paper we consider a simple example of quantum frustrated antiferromagnetism, 
namely the frustrated spin-1/2 Heisenberg model, with Hamiltonian 

(id) («J'> 

The spin operators obey Si ■ Si = 3/4, and Ji = 1 throughout this paper. The notations 
and indicate summation over the nearest- and next-nearest neighbor bonds 

on a square lattice, each bond being counted once. While the model has attracted most 
attention simplified mo del0 of the effects of doping on copper oxide planes in the 
high-temperature superconducting copper oxides, it is of rather more general interest. A 
complete understanding would provide a clear example of answers to several general ques- 
tions about quantum phase transitions. 

The first question is that even in a ground state with rather classical looking symmetry, 
in this case an antiferromagnet, how do we show unequivocally that the order really is of 
long range and not simply local? How do we calculate physically measurable correlations 
without relying on low order perturbation theory? In the present case, for small frustration 
the appearance, in the limit of infinite size, of spontaneous symmetry breaking is displayed 
in a relatively simple model. Indeed the renewed interest in the model was because of 
doubts that the unfrustrated case would display long-range order in the thermodynamic 
limit. While such doubts are now relatively rare thanks to extensive numerical calculations 
and tighter rigorous limits for higher spin and lower spin symmetry,! there is as yet no 
rigorous proof for the isotropic spin one-half model in two spatial dimensions. One reason 
for the present study is to test the quantitative success of ideas of finite size scaling as 
applied to numerical diagonalizations that are perforce limited to what seem unhelpfully 
small samples. 

The history of finite size effects goes back to Anderson in the nineteen-fifties,i who first 
invoked the fact that the infinite degeneracy of the ground-state with spontaneously broken 
continuous symmetry must be manifest in a large number of nearly degenerate states in a 
large but finite system. This idea of a "tower" of states whose degeneracy corresponds to 
the ultimate symmetry, and whose energy scales determine the long distance parameters of 



the spontaneously broken model of the infinite system has since been made more precise 
and less dependent on perturbative concepts in the language of non-linear sigma models.@ 
The model we consider here has the advantage over, for example, the triangular or Kagome 
antiferromagnetsi'i in that the classical limit has a simpler unit cell and thus the structure of 
the towers should be simpler to test. One of our aims here will be to show that it is possible 
to extract the parameters of the long wavelength physics in the ordered regime. In practice 
the difficulties of applying finite size studies are still considerable: there are subleading as 
well as leading corrections which make the ultimate goal of reliable quantitative calculations 
difficult even here. It is helpful that we may easily stabilize the ordered state to study the 
disappearance of order in a controlled fashion by applying negative J2. 

A second general question relevant to other quantum phase transitions, is whether the 
finite size methods developed can be applied all the way to a critical point at which the 
order may disappear with a continuous transition. The first step is to identify the parameter 
J2C of this critical point unequivocally; even its existence is still a matter for contention. 
Indeed some self-consistent spin-wave expansions have been interpreted as indicating a first 
order transitioni^, at least for large spin. We shall present results which we feel are rather 
convincing as to the existence of a critical point and a reasonably accurate estimate of its 
value. 

A third question, separate from the study of ordered antiferromagnetism, is the question 
of what happens when this order disappears. In the mapping of quantum interacting ground 
states to thermodynamics of classical models in higher dimension, there is at first sight a 
difference in that quantum phase transitions tend to show order-order rather than order- 
disorder transitions. Of course what one means by "order" is crucial to such a distinction. 
Here an ordered state would be understood to have long range order in a different local 
order parameter, for example a spin-Peierls dimerization variable or chirality parameter. 
In this paper we do not discuss in detail the nature of the intermediate state, but we do 
produce evidence that at least it corresponds to one of zero uniform susceptibility. This is 
compatible with either the chiral or dimerized phases favored previously.0 With increasing 
sample size both phases would have an exponentially vanishing ferromagnetic susceptibility. 



We use exact finite-size diagonalization on clusters of = 16, 20, 32, 36 sites (fig.|T]) with 
periodic boundary conditions. Tliese are all the clusters accessible by our present calcula- 
tional means that both respect the square symmetry of the lattice and do not frustrate the 
collinear magnetic state expected at large J2 (this last condition is violated, for example, 
for the 18 and 26 site clusters, and more generally whenever is not an integer multiple of 
4). The present study is an improvement of our previous finite-size calculation,^ which was 
restricted to the untilted four by four and six by six lattices. We shall see that the inclusion 
of the intermediate sizes is extremely useful in allowing us to test the consistency of the 
calculation. The final qualitative conclusions are not drastically changed but estimates of 
critical parameters in particular are altered. We are also able to articulate questions about 
whether it is advisable to extrapolate from the sixteen site cluster. In our previous study 
we had no choice even though the rather special hypercubic property of the 16 site cluster 
was a concern. 

As long as frustration is not too strong, the ground state of the model ( |1.1[ ) is expected 
to have long-range antiferromagnetic order, and then the low-energy long-wavelength exci- 
tations are expected to be described by the quantum nonlinear sigma model0, with action 

Here n = n(r,t) is the local orientation of the staggered magnetization, with |n| = 1, c 
and ps are the spin-wave velocity and spin stiffness, and the inverse temperature j3 has 
to be taken to infinity here as we are interested in ground state properties. Lowest order 



;i.2) 



spin- wave theory giveslllOco = a/2(1 — 2 J2/J1) Ji, pso = (-^1 — 2J2)/4, but there are of 
course important quantum fluctuation corrections to these quantities. One way to extract 



these corrections from finite size data will be discussed in sec. pi B| below. We note that 
the magnetic susceptibility at g = is given by x = Ps/(?i which in spin-wave theory 
equals 1/(8 Ji). A major aim of the paper is to obtain controlled estimates of the different 
parameters beyond spin-wave theory. A summary of current results for the unfrustrated 
case can be found in recent review articles.lii'lii 



II. NUMERICAL PROCEDURES AND RESULTS 



We wish to find eigenvalues and eigenvectors of tlie Hamiltonian ( |1.1|) on large clusters. 
In order to achieve this, and given that computational power is and will remain limited, it 
is necessary to use the symmetries of the problem to reduce the size of the corresponding 
Hilbert space as much as possible. For the = 16, 32, 36 clusters we use: 

1. translational symmetry (A^ operations for an A^-site cluster). 

2. reflection on horizontal (-R-) and vertical {R\) axes (4 operations). For the N = 16 
and = 36 cluster, both symmetry axis pass in between rows of spins. However, for 
A^ = 32, the i?_-axis coincides with the central row of spin (see fig.|l]). 

3. if a given eigenstate has the same eigenvalue under R_ and R\, then reflection on the 
diagonal running from the lower left to the upper right of the cluster (R/) is also a 
symmetry operation, and can be used to further reduce the size of the Hilbert space 
by a factor 2. For the 32 site cluster, this operation has to be followed by a translation 
to remap the cluster onto itself. 

4. if the z-component Sz of the total magnetization (which commutes with the Hamil- 
tonian) is zero, then the spin inversion operation | f) ^ | J.) is also a symmetry and 
leads to a further reduction by a factor 2. In principle a further considerable reduc- 
tion of the Hilbert space could be achieved by using the conservation of the total spin 
S"^. However, there does not seem to be any simple way to efficiently incorporate this 
symmetry. Note that we could equally have chosen reflection operators centered on a 
particular site. 

The point group operations Id, R_, R\, R/ generate the point group symmetry C^^. These 
operations are only compatible with the translational symmetry for states of momentum 
Q = or Q = (tt, 7r). In particular, for our clusters the ground state is always at Q = 0. For 
the 20 site clusters reflections are not symmetry operations, and we use rather a rotation 
by it/ 2 as generator of the point group. The symmetry group at the interesting momenta 
Q = or Q = (vr, 7r) then is C4. 



We use a basis set characterized by the value of Szi at each lattice site i. An up (down) 
spin is represented by a bit 1 (0) in a computer word. Thus, a typical spin configuration 
(e.g. for a linear system of 4 spins) would be represented as 

I TTiT) = 110I2 = 13 (2.1) 

To implement the symmetry, we do not work in this basis, but use rather symmetry-adapted 



basis states. E.g. to remain in the one-dimensional toy example, instead of (pTTI) we use 
the normalized basis state 

^(1 TTiT) + I TTTi) + I iTTT) + I TiTT)) = ^(|13) + |14) + |7) + |ii)) ^ |7) (2.2) 



where the lowest ( "minimal") integer of the 4 states occurring in ( |2.2| ) is used to represent 
the state. 

Our procedure to determine eigenvectors and eigenvalues procedes in three steps: (i) 
starting from an arbitrary basis state of given symmetry and S^, the whole Hilbert space 
is generated by repeated application of the Ji part of the Hamiltonian, and the basis set 
is stored; (ii) the Hamiltonian matrix is calculated and stored in two pieces, corresponding 
to the Ji and J2 parts of the Hamiltonian; (iii) the matrix is used in a Lanczos algorithm 
to obtain eigenvalues and eigenvectors of the Hamiltonian. The principal difficulty in steps 
(i) and (ii) is that application of the Hamiltonian to a state represented by a "minimal" 
integer will of course in general not produce another minimal integer state, but rather 
a state that needs to be brought into minimal form by the application of a symmetry 
operation. The trivial solution would be to try out all possible operations. This however 
would be extremely time consuming (there are 576 symmetry operations for the 36 site 
cluster!). Instead we use a different procedureEl: the basis states are coded in a computer 
word so that the R\ operation corresponds to the exchange of the two halfwords. Each 
halfword then can be an integer between and 2^/^. We then create a list specifying for 
each halfword the corresponding minimal state (integer) and the symmetry operation(s) 
connecting them. The length of this list is relatively moderate (2^^ at worst), and it can be 
easily kept in computer memory. 

The minimal state corresponding to a given basis state is now determined by looking in 
this list for the minimal states corresponding to the two halfwords. If necessary, the two 



halfwords are exchanged (i.e. a R\ operation is performed), so that the smallest of the two 
halfwords constitutes the high-bit halfword of the resulting state. Finally, the symmetry 
operation leading to this high-bit halfword is applied to the remaining halfword. In about 
80% of the cases this symmetry operation is uniquely determined. In the remaining cases, 
more than one symmetry operation has to be tried out in order to find the minimal state. 
However, the extra calculational effort is relatively small: e.g. for = 36, only for about a 
thousand out of 2^*^ possibilities are there more than eight symmetry operations to be tried 
out. Using this method for = 36 the CPU time needed to calculate the basis set and the 
Hamiltonian matrix is approximately 30 min and 90 min., respectively.^ For the smaller 
clusters, CPU time requirements are obviously much less. In table I we show the size of 
the basis set and the number of non-zero matrix elements for states of Ai {N = 16, 32, 36) 
or A {N = 20) symmetry at Q = 0. These are the subspaces containing the groundstate, 
apart from the case of relatively large J2 on the = 20, 36 clusters, where the ground 
state has point group symmetry B [N = 20) or Bi [N = 36). Note that the number of 



basis states for the larger clusters is very close to the naive expectation 



^ 1.57554 X 10^ for = 36). 

The number of matrix elements e.g. for A^ = 36 is still enormous. It is however obvious 
that the matrix is extremely sparse: on the average, there are fewer than 80 nonzero 
elements per line, which has in all 15, 804, 956 positions. One obviously only wants to 
store the addresses and values of the nonzero matrix elements. This would still need two 
computer words per non-zero matrix element, however, this requirement can be further 
reduced noting that all matrix elements are of the form Hij = Ji^2(Aj/Aj)/ij, where the 
Aj are the normalization factors of the symmetrized basis states (like the factor 1/2 in 
( |2.2D ), and the Jjj are small integers, which in the vast majority of cases equal unity. More 
specifically, Jjj is the number of times the action of the Hamiltonian on an unsymmetrized 
basis state \i) creates another unsymmetrized basis state \j) or a state related to \j) by a 
symmetry operation. The values of the Aj intervening in a given matrix element can be 
easily determined during the calculation, and we thus just store the positions of the unit 
integer. For the cases where /jj = n 7^ 1, the corresponding position is stored n times. 



Finally, a Cray computer word has 64 bits, and therefore can accomodate two addresses. 
In this way the whole matrix for the 36 site cluster can be stored in about 5 gigabytes, 
which is relatively easily available as disk space at the computing facility we are using. Space 
requirements could be further reduced by a factor 2 using the symmetry of the Hamiltonian 
matrix, however, this would have lead to a rather important loss of speed in the subsequent 
matrix diagonalization. 

To obtain the groundstate eigenvalue and eigenvector of the Hamiltonian we use the 
standard Lanczos algorithm, implemented by the Harwell library routine EA15AD. This 
routine performs rather extensive convergence checks and we thus avoid to perform unnec- 
essary time-consuming Lanczos iterations. The main problem at this level is the use of 
the still rather large matrix 5 gigabytes). The matrix clearly does not fit into the main 
memory of a Cray-2 (2 gigabytes). We therefore store the matrix on disk, and read it in 
by relatively small pieces, whenever a new piece is needed. This operation can be made 
computationally efficient by using "asynchronous" input operations, which allow one to 
perform calculations in parallel with the read-in operation for the next piece of the matrix. 
Moreover, using more than one input channel simultaneously the read-in operation can be 
further accelerated. In this way the total time overhang due to the continuous read-in of the 
matrix can be kept below 20% of total CPU time. To reach a relative accuracy of 10^^ for 
= 36, we need between 40 min. (J2 = 0) and 3 hours {J2/J1 ~ 0.6, slowest convergence) 
CPU time. We have performed a number of checks to insure the correctness of the numerical 
algorithm. The most important one is to calculate the groundstate energy for ferromagnetic 
interaction (Ji 2 < 0), which of course is known to be {Ji + J2)N/2. However, the numerical 
calculation in the Sz = subspace is nontrivial because the Hilbert space and, up to an 
overall minus sign, the matrix are of course the same as for the antiferromagnetic case. We 
also compared our results with previous finite size calculations0^0 (for = 16, 20, 32) and 
quantum Monte Carlo result (for A^ = 36, J2 = 0), and found agreement in all cases. Fi- 
nally, an independent check of the numerical accuracy of the Lanczos algorithm is provided 
by starting the Lanczos iterations with different initial vectors. In each case we found a 
relative accuracy of at least 10~^ for the ground state eigenvalues. Similarly, expectation 
values calculated with the eigenvector are found to have relative accuracy of 10^^. In table 



Ill we list ground state energies of the different clusters for a number of values of J2. A 
more complete set of results is displayed in fig.^. 

More important for the following analysis are the values of the Q-dependent magnetic 
susceptibility (or squared order parameter) 

MUQ) = ^^^^^^ 5^(5. ■ 5,)e"«-(^-^^) . (2.3) 



Following arguments by Bernu et al.i we use a normalization by a prefactor j^^^_^_2-) instead 
of the usual l/N"^. In the thermodynamic limit, these possibilities are obviously equivalent. 
However, for the relatively small cluster we are using, there are sizeable differences in the 
results of the finite-size scaling analysis. The choice in eq. ( [^.3[ ) is essentially motivated by 
the fact that in a perfect Neel state M^(Q) is entirely size-independent .0 More generally, 
this choice eliminates to a certain extent the overly strong contributions from the terms 
with i = j in eq. ( p.3|) . 

Some values of M^(Q) are shown in table |T|, and complete curves are in fig.§. The 
values displayed (and used in the following analysis) are always expectation values in the 
true ground state, e.g. for large J2 states of symmetry B [N = 20) or Bi {N = 36) are 
used. From the results shown it is quite obvious that the dominant type of magnetic order 
changes from Q = (vr, vr) at relatively small J2 (Neel state) to Q = (vr, 0) at larger J2 
(coUinear state). How exactly this change occurs will be clarified in the following section. 



III. FINITE SIZE SCALING ANALYSIS 

A. Order parameters 

The results shown in fig.|^ show a transition between a Neel ordered region for J2 < 0.5 to 
a state with so-called collinear order (i.e. ordering wavevector Q = (vr,0)) at J2 > 0.6. To 
analyze the way this transition occurs in more detail, we use finite-size scaling arguments.ii 
In particular, it is by now well established that the low-energy excitations in a Neel ordered 
state are well described by the nonlinear sigma model. From this one can then derive the 
finite-size properties of various physical quantities. The quantity of primary interest here 



is the staggered magnetization mo(Qo) defined by 

mo(Qo) = 2 hm Mn{Qq) , (3.1) 

where Qq = (vr, tt). The normahzation is chosen so that ?7io(Qo) = 1 in a perfect Neel state. 
The leading finite size corrections to niQ are given by 

MUQo) = ^mo(Qo)' + 1.2416^ + ... 



1 ,^ ,2/. 0.6208c , ,^ ^, 

4 PsVN 

where for the infinite system ki gives the amplitude of the diverging matrix element of the 
spin operator between the ground state and single magnon states at Q ~ Qo- 

Least square fits of our finite-size results to eq. ( |3.2| ) are shown in fig.^. For small 
values of J2 the scaling law is quite well satisfied: e.g. for J2 = the four data points in 
fig.|| very nearly lie on the ideal straight line, and the extrapolated value of the staggered 
magnetization, mo(Qo) = 0.649, is quite close to the best current estimates, mo(Qo) = 
0.615.0 Using the same type of finite size extrapolations for other values of J2, we obtain 
the results indicated by a dashed line in fig.^ . 

For J2 = 0, a check on the reliability of our method can be obtained by comparing 
the numerical results with what one would expect from eq.( p.2|) , using the rather reliable 
results for mo, c, and ps obtained by series expansion techniques.!!^ The curve expected 
from eq. ( |3.2|) is shown as a dash-dotted line in fig.§. It appears that there are sizeable but 
not prohibitively large next-to-leading corrections. 

Another measure of the reliability of the finite-size extrapolation can be obtained com- 
paring results obtained by the use of different groups of clusters. For negative J2, i.e. non- 
frustrating interaction, the values of mo(Qo) are nearly independent of the clusters sizes 
used, and the results in fig.|^ therefore are expected to be quite accurate. In this region 
the next nearest neighbor interaction stabilizes the antiferromagnetic order and therefore 
the staggered magnetization tends to its saturation value unity for large negative J2. On 
the other hand, for positive J2 the interaction is frustrating. In this case, the agreement 
between different extrapolations is less good. We note however, that in all but two cases 
the staggered magnetization tends to zero as in a second order phase transition, with a 



critical value of J2 between 0.34 and 0.6. The question than arises as to which extrapola- 
tion to trust most. In fact, none of the clusters considered here is free of some peculiarity: 
for = 16 and J2 = 0, there is an extra symmetry, because with only nearest neighbor 
interactions this cluster is in fact equivalent toa2x2x2x2 cluster on a four-dimensional 
hypercubic lattice; the = 20 cluster has a lower symmetry than all the others (C4 instead 
of Civ); for = 20 and = 36 the ground state changes symmetry with increasing J2; 
finally the 20 and 32 site clusters are unusual in that they are rotated, by different angles, 
with respect to the lattice directions. A priori, one might then argue that the best choice 
should be the least biased one, including all available clusters. As indicated by the dashed 
line in fig.|^, this leads to a critical value of J2 for the disappearance of antiferromagnetic 
order of ~ 0.48. 

However, from fig.^ it is quite clear that for J2 > 0.35 the 16 site cluster is highly 
anomalous in that M^(Qo) increases going to the next bigger cluster, whereas in all other 
cases there is a decrease with increasing size. Clearly, in fig.|^ a much better fit is obtained 
in this region by omitting the A^ = 16 results, leading to a reduced value, ~ 0.34 as 
indicated by the full line in fig.|^. The anomalous results obtained from the A^ = 16, 20, 32, 
A^ = 16, 32, 36 and A^ = 16, 32 fits are certainly due to an over-emphasis put onto the 
A^ = 16 results. Similar anomalous behavior of the 16 site cluster occurs in many cases 
in the region 0.3 < J2 < 0.8, and we therefore consider the results obtained using only 
A^ = 20, 32, 36 as more reliable. In particular, in this way we find a staggered magnetization 
of 0.622 at J2 = 0, only about one percent higher than the best current estimate, 'mo(Qo) = 
0.615. Beyond the precise value of the critical value J2C at which antiferromagnetic order 
disappears, the important result here, obtained by the majority of fits, is the existence of 
a second order transition, located in the interval 0.3 < J2 0.5. 

One might of course argue that it is not the A^ = 16 but rather the A^ = 20 cluster that 
is anomalous. However, closer inspection of the data in figs.^ and ^ clearly shows that the 
A^ = 20, 32, 36 data points remain reasonably well aligned even in the intermediate region 
0.3 < J2 < 0.8, whereas the alignment for A^ = 16, 32, 36 is much worse. In the following, 
we will therefore mostly rely on the A^ = 20, 32, 36 extrapolations. 



We now follow the same logic to analyze the behavior for larger J2, where fig.|^ suggest 
the existence of magnetic order with ordering wavevector Qi = (tt, 0). Of course, this 
state again breaks the continuous spin rotation invariance, and therefore the low energy 
excitations are described by a (possibly anisotropic) nonlinear sigma model. There is an 
additional breaking of the discrete lattice rotation symmetry (ordering wavevector (0, vr) is 
equally possible), however, this does not change the character of the low-lying excitations. 
The finite size behavior is entirely determined by the low energy properties, and therefore 
we expect a finite size formula analogous to eq.( p.2|) : 

MUQi) = lmo{Q^r + ^ + ... . (3.3) 



Here the factor 1/8 (instead of 1/4 in ( |3.2| )) is due to the extra discrete symmetry breaking 
which implies that finite-size ground states are linear combinations of a larger number of 
basis states. Moreover, the linear sigma model is anisotropic, because of the spontaneous 
discrete symmetry breaking of the ordering vector, and consequently a precise determination 
of the coefficient of the V^^term is not straightforward. The important point here is 
however the A^-dependence of the correction term in eq.( |3.3| ). 

Least square fits of our numerical results to eq. ( p.3| ) are shown in fig.||, and the extrap- 
olated collinear magnetization mo(Qi) is shown in fig.|^. For J2 > 0.8 eq.( ^73|) provides a 
satisfactory fit to our data, even though not quite as good as in the region J2 < in the 
staggered shown by the spread of different fits in fig.|^ (compare fig.^ in the region 

J2 < 0). For smaller J2 there is a wide spread in the extrapolated results, depending on 
the clusters used. We notice however that for the majority of clusters used, there is a com- 
mon feature: mo(Qi) remains finite down to J2 = 0.65, and then suddenly drops to zero 
at J2 = 0.6. This would indicate a first order transition to the collinear state somewhere 
in the interval 0.6 < < 0.65. This interpretation also seems consistent with the raw 
data fo fig.0: the increase of M^(Qi) around J2 = 0.6 is much steeper than the growth of 
M^{Qo) with decreasing J2. From the = 16, 20, 32, 36 extrapolation one then obtains a 
collinear magnetization which is roughly constant above J2C at mo((5i) ~ 0.6. Notice that 
the first-order character of the transition is not due to the level crossings occurring in the 
= 20 and iV = 36 clusters: if these clusters are omitted from the extrapolation, the first 



order character is in fact strongest (cf. fig.^. 

On the other hand, inclusion or not of the = 16 cluster plays an important role 
because this cluster shows again anomalous behavior in the region of intermediate J2: for 
J 2 < 0.7 Mff{Qi) increases when N increases from 16 to 20, contrary to what eg. ( p. 3D 
suggests. If one therefore omits the = 16 cluster from the extrapolation, results quite 
consistent with a second order transition are obtained, however this time with a critical 
coupling ~ 0.68. 

Beyond quantitative results, the most important conclusion of this analysis is the ex- 
istence of a finite interval without magnetic long range order: if all available clusters are 
included in the analysis, this interval is 0.48 ^ J2 ^ 0.6. if, because of the anoma- 
lies discussed above one omits the = 16 cluster, the nonmagnetic interval is increased 
to 0.34 < t/2 ^ 0.68. The study of the ground state symmetry in this region requires 
a detailed analysis of a number of different non-magnetic order parameters and will be 
reported in a subsequent paper. However, at this stage, the magnetic structure factor 
S{Q) = (A^ + 2)Mjf{Q) already gives some valuable information: in fact, as shown in fig.^, 
with increasing J2 the collinear peak at the X point grows and the Neel peak at the M point 
shrinks, however there never is a maximum at other points. There is thus no evidence for 
incommensurate magnetic order. 

B. Ground state energy, spin wave velocity, and stiffness constant 

The ground state energy per site in the thermodynamic limit can be obtained from the 
finite-size formula for an antiferromagnetS 

Eo(Ar)/iV = eo- 1.4372^ + ... , (3.4) 

where c is the spin-wave velocity. Again, in the collinear state, an analogous formula holds, 
but with c replaced by some anisotropy-averaged value. Fits of our numerical results are 
shown in fig.^ Away from the "critical" intermediate region, i.e. for J2 < 0.2 and J2 > 0.8, 
eq.( |3.4| ) provides a rather satisfying description of the results, in particular if the A^ = 16 
cluster is disregarded. The fit is even considerably better than that for the order parameters 
(compare fig.^. This is certainly in large part due to the much weaker finite size correction 



to the ground state energy, as compared to those for the order parameters. On the other 
hand, in the intermediate region 0.4 < J2 ^ 0.7, the fits are not very good. In this region the 
ground state energy per site is rather irregular, for example there is generally a decrease 
from = 32 to = 36 contrary to what eg. ( p.4[ ) suggests. The failure of eq.( p.4| ) in 
the intermediate region is of course not surprising, as the analysis of the previous section 
showed the absence of magnetic order, which implies the non-existence of an effective 
nonlinear sigma model and therefore the invalidity of the extrapolation formula (|3.4|) . The 
result of our extrapolations is shown in fig.|l^. Over most of the region shown, results 
from extrapolations using different clusters are indistinguishable on the scale of the figure. 
Only close to the critical region is there a spread of about 2 percent in the results. In 
particular, at J2 = we find values between cq = —0.668 and Cq = —0.670, very close to 
the probably best currently available estimate, obtained from large-scale quantum Monte 
Carlo calculations, of eo = —0.66934.11 



The amplitude of the leading correction term in eq.( |3.4| ) allows for a determination of 
the spin-wave velocity c. Results are shown in fig.|rT[ In this case, there is a wider spread 
in results. This is certainly not surprising, given that this quantity is derived from the 
correction term in eq.( |3.4| ). Nevertheless, the agreement between different extrapolations 
is reasonable for J2 < 0. At J2 = and using all clusters we find c = 1.44Ji, close to but 
somewhat lower than the best spin-wave result csw = l-65Ji. A smaller value is found 
from the A^ = 20, 32, 36 extrapolation: c = 1.28. For positive J2 the extrapolations give 
different answers, according to whether the A^ = 16 cluster is included or not. This of 
course is due to the anomalous behavior of this cluster in the energy extrapolations (see 
fig.^. An important point should however be noticed: independently of the inclusion of 
the A^ = 16 cluster, at the critical value J2C for the disappearance of the antiferromagnetic 
order the spin-wave velocity remains finite. 

The final parameter in the nonlinear sigma model is the spin stiffness constant ps- It 
can be found from our finite size resultsil 
rrioiQofc 

Ps = , (3.5) 



with Ki determined from eq.( |3.2| ) SB This relation determines the second form of eq.( p.2|) 
above. Results are shown in fig.|T2|. Again, for the same reasons as before, there is some 



scatter in the results, because of the use of the correction terms in eqs. ( p.2| ) and ( |3.4| ) . The 
results at J2 = {ps = 0.165 or 0.125 according to whether = 16 is included or not) is 
lower than other estimates (p^ ~ O.18Ji).0§i The fact that — as J2 ^ J2C is again in 
agreement with expectations from the nonlinear sigma model analysis, but is of course a 
trivial consequence of eq. (|3.5| ). 

In the collinear region, there is an additional anisotropy parameter in the effective 
nonlinear sigma model, and the corresponding effective parameters therefore cannot been 
obtained straightforwardly from the lowest finite size correction terms. 

C. Susceptibility 

An independent test of the reliability of our results can be obtained by calculating the 
susceptibility x'- even in an antiferromagnetically ordered state, the ferromagnetic suscep- 
tibility is finite, whereas for unconventional states (e.g. dimer or chiral), one has a spin 
gap and therefore a vanishing susceptibility. The vanishing of the susceptibility can thus 
be associated with the vanishing of the magnetic order parameter. Moreover, in an an- 
tiferromagnetic state one has x = Ps/c^i and we thus have a consistency check on our 
calculated values for c and p^. At fixed cluster size one has x{^) = where A^^ 

is the excitation energy of the lowest triplet state (which has momentum Q = (tt, tt) in 
an antiferromagnetic state). An extrapolation of x{^) to the thermodynamic limit can be 
performed using the finite-size formulaEllii x = x{^) ~ const./ y/N, and results are shown 
in fig.|l^. Again, the = 16 cluster behaves anomalously in that x(A^) increases going from 
A^ = 16 to A^ = 20, whereas for bigger clusters there is the expected decrease. In the present 
case, this anomaly occurs for nearly the whole range J2 > 0. Also, our result for J2 = and 
using A^ = 20,32,36 is x = 0.0671, very close to both Monte Carlo estimateJ^ and series 
expansion results.ii^ We therefore think that the A^ = 20, 32, 36 extrapolation is the most 
reliable one. It is rather pleasing to note that this independent estimate gives a critical 
value for the vanishing of the susceptibility (which indicates the disappearance of gapless 



magnetic excitations and therefore of long-range antiferromagnetic order) of J2c ~ 0.42, 
quite close to the estimate we found above by considering the order parameter. 

A quantitative comparison of results for the susceptibility obtained either from the 
excitation gap or from the previously calculated values of c and ps and using x = Ps/(? 
reveals considerable discrepancies (see fig.lTSI), even well away from the "critical region" 
J2 ~ 0.4. The most likely explanation for this is that our calculation of c and ps is based 
on corrections to the leading finite-size behavior, whereas x is obtained directly from the 
gap. In particular, judging from the case J2 = 0, we probably underestimate the spin wave 
velocity by quite a bit. The direct estimate of x is thus expected to be more precise. 

An analogous calculation of the susceptibility can be performed in the region of larger 
J2, where the lowest excited triplet state is at Q = (vr, 0). In this case, because of the 
double degeneracy of this state, the susceptibility is given by x = 2/(A^At). Because of 
the lower symmetry of the wavevector, the Hilbert space needed to determine the excited 
state roughly double in size, and for = 36 has dimension 31561400. We use the same 
finite-size extrapolation as before, and results obtained for different combinations of cluster 
sizes are shown in fig.|T^. The 16 site cluster again shows rather anomalous behavior and 
therefore we do not take it into account in these extrapolations. The results then indicate a 
transition into a nonmagnetic {x = 0) state at J2/J1 ^ 0.6, in approximate agreement with 
what we obtained from estimates of the order parameter above. The decrease of x with 
increasing J2 is not surprising, as for large J2 the model consists of two nearly decoupled 
unfrustrated but interpenetrating Heisenberg models, each with exchange constant J2, and 
consequently one has x l/<^2- What is a bit more surprising is the sharpness of the 
maximum of x around J2/J1 = 0.7. 



D. Comparison with spin— wave theory 

Linear spin-wave theory (LSWT) has proven to be a surprisingly accurate description 
of the ordered state of quantum antiferromagnets even for spin one-half. We here compare 
our numerical results with that approach. The lowest order spin-wave energies in the 



antiferromagnetic and collinear state are 

cOAFik) = 2{[1 - a{l - rj,)f - ^lY^' , (3.6) 
c^coiiik) = {{2a + 7,)2 - {2ar]k + lyfV^' , (3.7) 

where a = J2/J1, 7q = cos ka, 7fc = {ix + 7y)/2, and t]/. = 'jxly In LSWT, the order 
antiferromagnetic and colhnear order parameters then are given 

mo Qo = 1 - / d'k ^— 3.8 

moiQi) = 1 - A / d'k^-^ (3-9) 

where the integration is over the full first Brillouin zone. A comparison of our results with 
this approach is shown in fig.|T^. For the antiferromagnetic order parameter, we observe 
very satisfying agreement. What is slightly disturbing here is that inclusion of the next 
order {1/ S'^) correction to eq. actually makes the agreement worse, even for negative 

J2 where the next -nearest neighbor interaction stabilizes the order and spin-wave theory 
therefore should be increasingly reliable. For example, for J2 = —Ji these corrections lead 
to0 '"^o(Qo) = 0.775, whereas we find mo(Qo) = 0.846. To which extent higher-order spin 
wave theory can be used systematically even in this region thus seems unclear. For the 
more interesting case of positive J2, higher corrections to spin-wave theory give more and 
more strongly diverging results as J2 — * 0.5, and it is not clear how any useful information 
can be obtained from these higher order corrections in the frustrated case. We therefore 
limit our comparison here to linear spin-wave theory. 

For the collinear state at large J2, there is a similar good agreement between spin- 
wave theory and our results, except for the immediate vicinity of the transition to the 
nonmagnetic state. Moreover, it appears that the prder parameter tends for large J2 to a 
value very close or identical to that of the antiferromagnetic order parameter at J2 = 0. This 
is in fact not difficult to understand: for J2 ^ Ji our model represents two very weakly 
coupled sublattices, with a strong antiferromagnetic coupling J2 within each sublattice. 
Consequently, the ground state wave function is to lowest order in J1/J2 a product of the 
wavefunctions of unfrustrated Heisenberg antiferromagnets on the two sublattices. We then 
obtain M^(Qi, J2 = 00) = Mj^^^{Qo, J2 = 0)/2, and thus from eqs. (|3.2|) and ( p.3|) we have 



the exact result 



hm mo(Qi) = mo(Qo)|j2=o • (3.10) 

J2/Jl—>00 



The spin-wave results ( |3.(j| ) and ( |3.9| ) as well as our numerical results do satisfy this relation. 
To obtain this agreement it was however crucial to use the normalization of Mj^{Q) shown 
in eq.( |2.3| ). Using a factor l/N"^ instead we obtain for large J2 ^o(Qi) ~ 0.4, which is far 
too low. The reason for this is that our extrapolation with = 16, 20, 32, 36 corresponds, 
for large J2, to a calculation on two nearly uncoupled and unfrustrated sublattices, each 
with = 8, 10, 16, 18. On such small lattices, short-range effects are obviously rather 
large, and therefore the proper normalization of is particularly important. 

The ground state energy per site is given in lowest order spin-wave theory by0'il 

eo = |(« - 1) + /" (fk iUAFik) , (3.11) 



3 1 

eo = + ^ I d''k uj,ouik) (3.12) 



2'" ' 87r2 
2" ' 87^ 

for the antiferromagnetic and collinear state, respectively. As can be seen in fig.|l3, these 
results are rather close to our finite-size extrapolations. Nevertheless, there is a significant 
discrepancy: e.g. for J2 = the spin-wave result is eo = —0.6579, compared to the 
presumably best estimate from large scale Monte Carlo calculationsil, Cq = —0.66934. On 
the other hand, as discussed above, our finite size extrapolation gives values very close to 
this. It would thus seem that, as far as the ground state energy is concerned, finite-size 
extrapolation is more precise than linear spin-wave theory. 



Comparing our results for the spin-wave velocity and spin stiffness (figures |Tl| and |12|) 
to the LSWT results c = ^2(1 - 2J2/J1) Ji and = (Ji -2J2)/4, one finds rather sizeable 
discrepancies, both for ps and for c. Nevertheless, the functional form for large negative J2 
seems to be correct. However, here a detailed comparison seems not particularly useful as 
LSWT results themselves are rather imprecise (as shown e.g. by the large renormalization 
of the susceptibility at J2 = 0). 



IV. SUMMARY AND DISCUSSION 



In this paper we have reported detailed finite-size calculations on the frustrated spin-1/2 
antiferromagnetic Heisenberg model on the square lattice. Using finite-size extrapolation 
formulae, we derived results for a number of physical properties. The most important find- 
ing seems to be the existence of a region of intermediate second nearest neighbor coupling 
J2 where no magnetic order, antiferromagnetic, collinear or otherwise, exists. The location 
of the boundaries of this nonmagnetic region depends on the cluster size involved in the 
estimate. For = 16,20,32,36 we find the interval 0.48 < J2/ Ji < 0.6 to be nonmag- 
netic, whereas with = 20,32,36 this interval is larger: 0.34 < J2/J1 < 0.68. Given the 
irregular behavior of the = 16 cluster we often found above, in particular in the region 
of intermediate J2, the second estimate would appear to be the more reliable one. In any 
case, independently of which extrapolation one prefers, there is a nonmagnetic interval. 

Beyond the existence of a nonmagnetic region, we have also obtained quantitative esti- 
mates for a number of fundamental physical parameters in the magnetically ordered states, 
antiferromagnetic for small or negative J2, collinear for large positive J2. The accuracy 
of these estimates can best be assessed by comparing with the unfrustrated case J2 = 0, 
for which case there are currently rather precise results available, mainly from large-scale 
Monte Carlo calculations and series expansions. A summary of our results, together with 
other recent data, is given in table |^ Our results for the ground state energy, the antifer- 
romagnetic order parameter, and the susceptibility agree to within a percent or better, with 
the best currently available numbers. Finally, our estimates for the spin-wave velocity and 
the spin stiffness are rather imprecise. This is certainly mainly due to the fact that these 
quantities are obtained from the amplitudes of the leading correction to the asymptotic 
large-size behavior of the ground state energy and the order parameter susceptibility, and 
these correction are almost certainly estimated less precisely than the leading terms. 

We found it instructive to also investigate regions where magnetic order is well- 
established, i.e. J2 ^ for the antiferromagnetic case and J2 > Ji for the collinear case. In 
these regions we find that the finite-size formulae like ( |3.2| ) and ( |3.4| ) provide an excellent 
fit to our numerical results. The progressive worsening of the quality of the fits as the in- 



termediate region is approached certainly is consistent with the existence of a quahtatively 
different ground state in that region. If on the contrary the transition between antiferro- 
magnetic and colhnear order occurred via a strong first order transition (as suggested by 
some approximate theories, see below), no such progressive worsening is expected. We also 
notice in this context that the = 16 cluster is systematically the one exhibiting the largest 
deviations from the expected behavior, probably due to its unusually high symmetry. We 
thus feel that estimates ignoring this cluster may be more reliable. 

Another way to assess the consistency of the finite-size extrapolations we are using is 
to verify the underlying scaling hypothesis via a "scaling plot" . The fundamental constants 
c and ps of the nonlinear sigma model define a length scale c/ ps-, and if finite size scaling 
is verified one therefore expects all finite size corrections to be universal functions of the 
variable x = c/{psVN). In particular, for the order parameter susceptibility we expect 

MUQo) = mo{Qof<!>{x) . (4.1) 



Combining the second form of eq.( |3.2| ), eq. ( |3.5| ), and this definition, the small-x expansion 
of the scaling function is $(x) = (1 + 0.6208x)/4. Plots of our results for M^(Qo) as a 
function of the scaling vaiable x are shown in fig.|16[ One sees that for the = 20, 32, 36 
extrapolation the plot is nearly perfect in that nearly all data points are collapsed onto 
a single curve. The only points that show a significant deviation are those obtained for 
iV = 16 close to the phase transition to the nonmagnetic state. This of course is nothing 
but a manifestation of the anomalous behavior of this cluster already found previously. The 
behavior for the = 16, 20, 32, 36 extrapolation is clearly less satisfying. A similar scaling 
plot for the ground state energy produces even better results, due to the better convergence 



of the corresponding finite-size formula ( |3.4| ). 



A scaling plot like fig.|T6| permits to assess the consistency of data obtained for clusters 
of different sizes, however, the form of the scaling function itself is obviously less significant 
as the coefficients c and ps entering the definition of the scaling variable x are calculated 
assuming finite-size formulae like (|3.2| ) and ( p.4|) to be valid, i.e. implicitly assuming the 
form $(x) = (1 + 0.6208x)/4. An independent estimate of $ can in principle be obtained 
using independent estimates for c and pg. We do not have currently such an estimate for 



Ps, however we can use our independent results for the susceptibihty (fig.0) to rewrite the 



scahng variable as x = 1/ y/xplN. The plot obtained using estimates for ps and x from 
N = 20, 32, 36 is shown in fig.p!7|. The collapse of data obtained for different sizes and 
values of J2 is not as satisfactory as in the previous case, however, this is certainly related 
to the fact that here we use a second independently estimated quantity, namely x- Still, for 
X < 5, the collapse is rather good, showing the consistency of our analysis in this region. 
For the larger clusters, this region corresponds to J2/J1 < 0.25, i.e. it extends rather 
close to the transition which occurs at J2/J1 ~ 0.34. For small x the calculated scaling 



function essentially agrees with the spin-wave results shown by the dashed line in fig.[T7. 
For X > 5, there are discrepancies between results obtained from M^(Qo) for different A^. 
This probably indicates that at least for the smaller clusters, finite size effects become so 
important that it is no more sufficient to include the lowest order finite size corrections 
only. The fact that the numerically found scaling function is larger than the spin-wave 
approximation is not entirely unexpected: in fact, for large x, i.e. in the critical region, one 
would expect ^(x) oc x^'^'^, where 77 is the correlation exponent of the three-dimensional 
Heisenberg model. However, we doubt that what we observe in fig.O is actually a critical 
effect. First, the numerical valueS of = 2 — (7///) is very small: r] 0.03, and one thus 
expects an extremely smooth crossover. Moreover, in fig.|l3 we have used the independently 
calculated susceptibility (see fig(|l3D which goes to zero only at J2/J1 ~ 0.42, rather than 
at J2/J1 ~ 0.34 where our estimated staggered magnetization vanishes. Consequently, the 



abscissae of the data points in fig.|T^ are underestimated, i.e. the data in fig.|T3 overestimate 
the true $(x). 

The Ji — J2 model we have studied here has been investigated previously by number 
of techniques. Previous finite-size studiesSIll found some indication of an intermediate 
phase without magnetic order, however due to the limitation to = 16 and 20 only, it was 
impossible to make extrapolations to the thermodynamic limit and to arrive at quantitative 
statements. Our own previous studyS, using = 16 and 36, produced results very similar 
to our current best estimates. However, due to the larger number of clusters we now use 
(and due to the possibility to ignore the anomalous = 16 cluster), we feel that our 
conclusions are considerably more reliable. 



Lowest order spin-wave theoryil produces a phase diagram very similar to ours (see 



fig.|l5]). On the other hand, higher order (in 1/5*) calculations do not seem to be very 
useful, due to increasingly strong singularities at J2 = -^i/2. It has been attempted to 
include higher order corrections using a selfconsistently modified spin-wave theoryJJi These 
calculations as well as the closely related Schwinger boson approachU produce a first order 
transition between Neel and coUinear state. A combination of Schwinger boson results and 
a renormalization group calculation^ gives on the other hand a second order transition 
from the Neel state to a magnetically disordered state, at J2c/ Ji = 0.15.0 However, the 
applicability of these approaches to an 5 = 1/2 system is hard to judge, mainly due to the 
absence of a small parameter that would make a systematic expansion possible. 

Quantum Monte Carlo methods are plagued with the sign problem for frustrated spin 
systems. Nevertheless, conclusions very similar to the modified spin wave calculations have 
been reached recently using a quantum Monte Carlo method.0 However, these results have 
rather large error bars and in some cases, in particular in the region of intermediate J2, are 
in disagreement with our present exact results. The validity of these results thus appears 
doubtful to us. 

Another approach has been via series expansion methods around a lattice covered by 
isolated dimers.0 Expanding around a columnar arrangement of dimers, these authors 
find a phase diagram very similar to ours, at least as far as magnetic order is concerned. 
However, these results are not without ambiguity: expanding around a staggered dimer 
arrangement, there appears to be a first order transition between Neel and coUinear states. 
The results of this method thus appear to be biased by the starting point of the expansion. 

The most obvious question left open by the present study is the nature of the ground 
state in the intermediate nonmagnetic region. Work extending our previous analysis^ is 
in progress and will be reported in a subsequent publication. It would also be interesting 
to investigate dynamical correlations functions, in particular in the vicinity of the critical 
point of the Neel state, = 0.34. One thus might gain additional insight into dynamical 
properties at a quantum critical point .00 Finally, one might also try to extend the size of 
the available clusters, in order to achieve better accuracy and reliability. The next useful 
cluster has 40 sites, and should be tractable in the near future. However, the next step then 



would be a cluster of 52 sites which would require computational means both in memory 
size and speed three or four orders of magnitude more powerful than what is currently 
available. A viable alternative to increase the size of the tractable clusters might be to 
combine the exact solution of moderately big clusters with Monte Carlo type approaches. 
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TABLES 

TABLE I. The number of states in the Hilbert space (uh) and the number of nonzero 

off-diagonal matrix elements (ne) for the clusters used in this paper. The numbers are for states 
in the Ai representation {A representation for N = 20) at momentum Q = 0. 



N Uh Tie 

16 107 3664 

20 1,321 55,660 

32 1,184,480 78,251,988 

36 15,804,956 1,170,496,152 



TABLE II. The ground state energy per site for different clusters and different values of J2 
(Ji is normalized to unity). Where the ground state representation changes with increasing J2 
(A^ = 20,36), the energies of both relevant representations are given. Boldface indicates the 
approximate location of changes in the ground state symmetry. 



J2 


16 


20(^) 


20(5) 


32 


36(^1 ) 


36(Bi) 


-1.00 


-1.16457 


-1.15103 


-0.801770 


-1.13251 


-1.12922 




-0.50 


-0.927249 


-0.915408 


-0.648275 


-0.900134 


-0.897626 




0.00 


-0.701780 


-0.690808 


-0.519508 


-0.680179 


-0.678872 


-0.603912 


0.10 


-0.659817 


-0.648444 


-0.501316 


-0.639048 


-0.638096 




0.20 


-0.619874 


-0.607519 


-0.487925 


-0.599542 


-0.599046 




0.30 


-0.582984 


-0.568545 


-0.479923 


-0.562283 


-0.562459 




0.40 


-0.551147 


-0.532381 


-0.476480 


-0.528379 


-0.529745 




0.50 


-0.528620 


-0.500615 


-0.476624 


-0.500096 


-0.503810 


-0.493941 


0.55 


-0.523594 


-0.487338 


-0.478122 


-0.489517 


-0.495178 


-0.490396 


0.60 


-0.525896 


-0.491633 


-0.491816 


-0.484599 


-0.493239 


-0.492267 


0.65 


-0.539382 


-0.516444 


-0.517029 


-0.502147 


-0.506588 


-0.506582 


0.70 


-0.563858 


-0.543309 


-0.545677 


-0.527741 


-0.529951 


-0.530001 


0.80 


-0.627335 


-0.600092 


-0.609595 


-0.586871 


-0.585428 


-0.586487 


0.90 


-0.696866 


-0.659162 


-0.677703 


-0.651509 


-0.645445 


-0.649052 


1.00 


-0.768468 


-0.719583 


-0.747576 


-0.718414 


-0.707495 


-0.714360 


1.20 


-0.914286 


-0.842827 


-0.88967 


-0.854910 




-0.848364 


1.50 


-1.13578 


-1.03098 


-1.10536 


-1.06229 




-1.05268 


2.00 


-1.50771 


-1.34863 


-1.46744 


-1.41044 




-1.39633 



TABLE III. The normalized susceptibility (eq.(|2 
clusters and different values of J2I J\- 



I)) at Q = (vr, vr) and Q = (vr, 0) for different 



J2/J: 



16 



M(7r,7r) 
20 32 



36 



16 



M(7r,0) 
20 32 



36 



-1.00 0.26924 0.26002 0.24296 0.23943 0.02778 

-0.50 0.26297 0.25221 0.23131 0.22660 0.02780 

0.00 0.24580 0.23430 0.20621 0.19879 0.02789 

0.10 0.23853 0.22785 0.19745 0.18893 0.02798 

0.20 0.22811 0.21949 0.18616 0.17601 0.02818 

0.30 0.21212 0.20816 0.17090 0.15800 0.02868 

0.40 0.18589 0.19193 0.14887 0.13109 0.03031 

0.50 0.14236 0.16693 0.11487 0.09236 0.03709 

0.55 0.11276 0.14834 0.09165 0.07062 0.04771 

0.60 0.07819 0.02915 0.05113 0.04378 0.07154 

0.65 0.04290 0.02015 0.01692 0.01954 0.10897 

0.70 0.02092 0.01303 0.01161 0.01232 0.13598 

0.80 0.00721 0.00561 0.00520 0.00611 0.15407 

0.90 0.00374 0.00302 0.00251 0.00314 0.15930 

1.00 0.00236 0.00193 0.00147 0.00183 0.16164 

1.20 0.00126 0.00103 0.00072 0.00088 0.16379 

1.50 0.00067 0.00055 0.00037 0.00044 0.16507 

2.00 0.00033 0.00027 0.00018 0.00021 0.16586 



0.02273 0.01470 0.01316 

0.02273 0.01471 0.01316 

0.02278 0.01476 0.01322 

0.02284 0.01480 0.01326 

0.02292 0.01489 0.01335 

0.02309 0.01506 0.01354 

0.02348 0.01545 0.01404 

0.02452 0.01669 0.01594 

0.02621 0.01880 0.01965 

0.11508 0.04627 0.03822 

0.12615 0.10333 0.08167 

0.13461 0.11321 0.10006 

0.14383 0.12265 0.11370 

0.14759 0.12616 0.11925 

0.14944 0.12759 0.12183 

0.15118 0.12876 0.12418 

0.15227 0.12939 0.12553 

0.15295 0.12977 0.12637 



TABLE IV. Comparison of our results at J2 = obtained from the A'^ = 16, 20, 32, 36 and 
N = 20, 32, 36 extrapolations with previous estimates from series expansions and quantum Monte 
Carlo calculations. A more complete compilation of previous results can be found in review 
articles 



eo 



mo{Qo) 



X 



N = 16,20,32,36 
N = 20,32,36 
series expansions ^ 
quantum Monte Carlo ^ 



-0.6688 
-0.6702 
-0.6696 
-0.6693 



0.649 
0.622 
0.614 
0.615 



0.0740 
0.0671 
0.0659 
0.0669 



^See refs. ^ and |27 
^See refs. |2l and ^ 



FIGURES 



FIG. 1. The clusters used in this paper. 



FIG. 2. The ground state energy per site as a function of J2/J1 for A'^ = 16 (full line), N = 20 
(dashed line), = 32 (dash-dotted line), and = 36 (dotted line). For clarity, the curves for 
N = 20, 32, 36 are also displayed shifted upwards by 0.2, 0.4, and 0.6, respectively. For = 16, 20 
we have results for J2/J1 in steps of 0.01, and only a continuous curve is displayed. For N = 32, 36, 
we have only results at the points indicated, and lines are a guide to the eye. 



FIG. 3. The magnetic susceptibility M{Q) at (a) Q = (vr, vr) and (b) Q = (vr, 0). The symbols 
and linetypes are the same as in fig.|2[ 

FIG. 4. Finite size results for M^(Qq) for different values of J2. The dashed lines are least 
squares fits to the data according to eq.(3.2), using all available clusters. The full lines are fits 
using only N = 20, 32, 36. The dash-dotted line is the leading finite size behavior expected at 
J2 = (see eq.(3.2)). 



FIG. 5. The staggered magnetization mo(Qo) ^ function of J2/J1 using different combina- 
tions of clusters (a). In (b) the "critical" region J2 > is shown enlarged. 

FIG. 6. Finite size results for M'^{Qi) for different values of J2. The dashed lines are least 
squares fits to the data according to eq.(3.3), using all available clusters. The full lines are fits 
using only TV = 20, 32, 36. 



FIG. 7. The collinear magnetization mo(Qi) as a function of J2/J1 using different combina- 
tions of clusters (a). In (b) the "critical" region 0.5 < J2 < 1-0 is shown enlarged. 

FIG. 8. Magnetic structure factor, as obtained from the = 36 cluster, in the Brillouin zone 
for J2/Ji=0 (•), 0.55 (A), 0.6 (o), 0.65 (★), 1 (■). The points F, M, X are Q = 0,Qo,Qi, 
respectively. Note that nowhere there is a maximum at a point different from M or X. 



FIG. 9. Finite size results for the ground state energy per site for different values of J2. The 
full lines are least squares fits to the data according to eq.(3.4), using all available clusters. The 
dashed lines are fits using only N = 20, 32, 36. 

FIG. 10. Ground state energy per site as obtained from finite size extrapolation using eq.(3.4). 
In the intermediate region 0.4 < J2 < 0.65 the extrapolation can not be used reliably, and no 
results are shown. Results obtained using different clusters are undistinguishable on the scale of 
this figure. The dash-dotted line is the spin-wave result, eqs.(3.11) and (3.12). 

FIG. 11. The spin wave velocity in the antiferromagnetic state as obtained from finite size 
extrapolation using eq.(3.4). No results are shown in the region where according to the previous 
analysis there is no antiferromagnetic order (J2 > 0.48 or J2 > 0.34 according to whether the 
iV = 16 cluster is included or not). 



FIG. 12. The spin stiffness in tlie antiferromagnetic state as obtained from finite size extrap- 
olation using eq.(3.5). Lines are a guide to tlie eye. 

FIG. 13. The susceptibility in the Neel region as obtained from x = 1/(-^At) (circles and 
triangles) and from x = Ps/c^ (crosses) using different extrapolations. As discussed in the text, 
the = 20, 32, 36 extrapolation is expected to be the most reliable one. 

FIG. 14. The susceptibility in the collinear region as obtained from x = 2/(AAj') using 
different extrapolations. As discussed in the text, the N = 20, 32, 36 extrapolation is expected to 
be the most reliable one. 

FIG. 15. Comparison of our finite size fits for the antiferromagnetic and collinear order pa- 
rameters (left and right curves, respectively) with linear spin wave theory. 

FIG. 16. Scaling plot of ^>(x) = MI{Qq) /mQiQ^f function of the variable x = c/ {ps^fN), 
using the N = 20, 32, 36 (lower curve) and N = 16, 20, 32, 36 (upper curve) extrapolations for c 
and ps- For clarity, data for the N = 16,20,32,36 extrapolation are shifted upward by 3 units. 
The straight lines represent the spin wave result ^{x) = (1 -|- 0.6208x)/4 

FIG. 17. Scaling plot of $(x) = MI{Qq) /mQ{QQf function of the variable 

X = 1/{xPsNY/'^, using the N = 20,32,36 results for x (see fig.||) and the N = 20,32,36 
extrapolation for p^. The dashed line represents the spin wave result ^{x) = (1 -|- 0.6208x)/4 
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N = 32 



N = 36 



FIG. 1. The clusters used in this paper. 
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FIG. 2. The ground state energy per site as a function of J2/J1 for A'" = 16 (full line), N = 20 
(dashed line), A'' = 32 (dash-dotted line), and N = 36 (dotted line). For clarity, the curves for 
N = 20, 32, 36 are also displayed shifted upwards by 0.2, 0.4, and 0.6, respectively. For N = 16, 20 
we have results for J2/J1 in steps of 0.01, and only a continuous curve is displayed. For N = 32, 36, 
we have only results at the points indicated, and lines are a guide to the eye. 
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FIG. 3. The magnetic susceptibility M{Q) at (a) Q = (vr, tt) and (b) Q = (vr, 0). The symbols 
and linetypes are the same as in fig. 2. 
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FIG. 4. Finite size results for M^(Qq) for different values of J2. The dashed Hues are least 
squares fits to the data according to eq.(3.2), using all available clusters. The full lines are fits 
using only = 20,32,36. The dash-dotted line is the leading finite size behavior expected at 
J2 = (see eq.(3.2)). 
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FIG. 5. The staggered magnetization mo((5o) function of J2/J1 using different combina- 
tions of clusters (a). In (b) the "critical" region J2 > is shown enlarged. 
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FIG. 6. Finite size results for M^{Qi) for different values of J2. The dashed Hues are least 
squares fits to the data according to eq.(3.3), using all available clusters. The full lines are fits 
using only N = 20,32,36. 



36 



(a) 




0.5 0.6 0.7 0.8 0.9 1.0 



h/h 

(b) 

FIG. 7. The collinear magnetization mo(Q^) as a function of J2/J1 using different combina- 
tions of clusters (a). In (b) the "critical" region 0.5 < J2 < 1.0 is shown enlarged. 
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FIG. 8. Magnetic structure factor, as obtained from the N = 36 cluster, in the Brillouin zone 
for J2/Ji=0 (•), 0.55 (A), 0.6 (o), 0.65 1 (■). The points T, M, X are Q = 0,Qo,Qi, 
respectively. Note that nowhere there is a maximum at a point different from M or X. 
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FIG. 9. Finite size results for the ground state energy per site for different values of J2. The 
full lines are least squares fits to the data according to eq.(3.4), using all available clusters. The 
dashed lines are fits using only N = 20, 32, 36. 
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FIG. 10. Ground state energy per site as obtained from finite size extrapolation using eq.(3.4). 
In the intermediate region 0.4 < J2 < 0.65 the extrapolation can not be used reliably, and no 
results are shown. Results obtained using different clusters are undistinguishable on the scale of 
this figure. The dash-dotted line is the spin-wave result, eqs.(3.11) and (3.12). 
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FIG. 11. The spin wave velocity in the antiferromagnetic state as obtained from finite size 
extrapolation using eq.(3.4). No results are shown in the region where according to the previous 
analysis there is no antiferromagnetic order (J2 > 0.48 or J2 > 0.34 according to whether the 
A/^ = 16 cluster is included or not). 
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FIG. 12. The spin stiffness in the antiferromagnetic state as obtained from finite size extrap- 
olation using eq.(3.5). Lines are a guide to the eye. 




FIG. 13. The susceptibiHty in the Neel region as obtained from x = l/(-^^r) (circles and 
triangles) and from x = Ps/c^ (crosses) using different extrapolations. As discussed in the text, 
the N = 20, 32, 36 extrapolation is expected to be the most reliable one. 
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FIG. 14. The susceptibility in the collinear region as obtained from x = 2/(A?^At) using 
different extrapolations. As discussed in the text, the N = 20,32,36 extrapolation is expected to 
be the most reliable one. 
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FIG. 15. Comparison of our finite size fits for the antiferromagnetic and collinear order pa- 
rameters (left and right curves, respectively) with linear spin wave theory. 
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FIG. 16. Scaling plot of = M|r(Qo)/mo(Qo)^ 

as a function of the variable x = c/{psVN), 
using the N = 20,32,36 (lower curve) and N = 16,20,32,36 (upper curve) extrapolations for c 
and ps- For clarity, data for the N = 16,20,32,36 extrapolation are shifted upward by 3 units. 
The straight lines represent the spin wave result ^{x) = (1 + 0.6208a;)/4 
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FIG. 17. Scaling plot of <^{x) = M|r(Qo)/"''o(Qo)^ ^ function of the variable 
X = l/(x/9siV)^/^ using the N = 20,32,36 results for x (see fig.l3) and the N = 20,32,36 
extrapolation for pg. The dashed line represents the spin wave result ^{x) = (1 + 0.6208a;)/4 
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